#Stuart's spatial poisson model, upscaled to populations, with no evolution.
# as implemented by Ben
#Last updated 10 November 09

# Build and load C functions
setwd("/Users/ben/Documents/Papers/Artificial waterbodies") #set the working directory


#####-------------------Functions-----------------#####
#The kernel
dcncross<-function(x, u, v) {  #stuart's cauchy-normal distribution in 2D
	(u^v*v*sqrt(v^v*(u^2*v+x^2)^(-2-v)))/(2*pi)
} 

#Loads spatial opportunities data with initial presences and
# 	site specific kernel parameters (u and v)
init.pop<-function(file.name){
	d<-read.table(file.name, header=T, sep="/t")
	names(d)<-c("ID","X", "Y", "Pres", "n.pairs", "u", "v")
	d
}

#loads pairwise distances from an R object
# object is a list of length (n points) with each list object being a dataframe with snk.ID and distance
init.pair<-function(file.name){
	d<-load(file.name)
	d
}

# spreads the population over gen generations
spread<-function(pop, gens, pairs, K){ #pairs is a list from init.pair
	for (i in 1:gens){
		plotter(pop)
		occp<-subset(pop, pop$Pres==1 & pop$age<6) #collect occupied sites
		potl<-pairs[occp$ID] #collect relevant parts of pair list
		potl<-do.call("rbind", potl)
		src.ID<-rep(occp$ID, times=occp$n.pairs)
		potl<-cbind(src.ID, potl)
		U<-rep(occp$u, times=occp$n.pairs) #expand source specific kernel parameters
		V<-rep(occp$v, times=occp$n.pairs)
		recruits<-dcncross(potl$dist, U, V) #calculate probabilities of arrival for each pair
		n.inds<-sum(recruits)*K #Calculates number of individuals to draw
		#browser()
		recruits<-tapply(recruits, potl$snk.ID, sum) #vector of probabilities for the multinomial
		recruit.ID<-as.integer(row.names(recruits))
		recruits<-rmultinom(1, n.inds, recruits)
		recruits<-data.frame(recruit.ID, recruits)
		recruits<-subset(recruits, recruits>2)
		pop$Pres[match(recruits$recruit.ID, pop$ID)]<-1
		pop$age[which(pop$Pres==1)]<-1+pop$age[which(pop$Pres==1)]
	}
	pop	
}

#plots individuals and a histogram of grwoth values in the population.
plotter<-function(popmatrix){
	plot(popmatrix[,2], popmatrix[,3], xlab="x-axis (arbitrary units)", ylab="y-axis (ditto)", pch=19)
	occp<-subset(popmatrix, popmatrix$Pres==1)
	points(occp[,2], occp[,3], pch=21, col="red")
}

#function for generating test data table of pairwise distances
pdist<-function(X, Y, maximum){
	out<-vector("list", length=length(X))
	for (i in 1:length(X)){
		xdiff<-(X-X[i])^2
		ydiff<-(Y-Y[i])^2
		dist<-sqrt(xdiff+ydiff)
		#src.ID<-rep(i, length(X))
		snk.ID<-1:length(X)
		part<-data.frame(snk.ID, dist)
		part<-subset(part, part[,2]<=maximum)	
		out[[i]]<-part
	}	
	out
}

##---------Test data-----------------##
#takes about 48 minutes to process 100,000 pairwise points on Ben's MacBook
n<-1000
X<-runif(n, 0, 1000)
Y<-runif(n, 0, 1000)
age<-rep(0, 1000)
ID<-1:n
Pres<-c(1, rep(0, n-1))
prs<-pdist(X, Y, 100)
n.pairs<-do.call("c", lapply(prs, nrow))
u<-rep(10, 1000)
v<-rep(3, 1000)
tabl<-as.data.frame(cbind(ID, X, Y, Pres, n.pairs, u, v, age))


spread(tabl, 10, prs, 100000)


######-------------------The Model-----------------#####

